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Identification  of  FIR  Wiener  systems  with  unknown,  non-invertible,  polynomial  non-linearities 

SETH  L.  LACYf*  and  DENNIS  S.  BERNSTEIN  J 

Wiener  systems  consist  of  a  linear  dynamic  system  whose  output  is  measured  through  a  static  non-linearity.  In  this  paper 
we  study  the  identification  of  single-input  single-output  Wiener  systems  with  finite  impulse  response  dynamics  and 
polynomial  output  non-linearities.  Using  multi-index  notation,  we  solve  a  least  squares  problem  to  estimate  products 
of  the  coefficients  of  the  non-linearity  and  the  impulse  response  of  the  linear  system.  We  then  consider  four  methods  for 
extracting  the  coefficients  of  the  non-linearity  and  impulse  response:  direct  algebraic  solution,  singular  value  decomposi¬ 
tion,  multi-dimensional  singular  value  decomposition  and  prediction  error  optimization. 


1.  Introduction 

Non-linear  system  identification  remains  one  of  the 
most  challenging  and  potentially  useful  problem  areas  in 
system  theory.  Numerous  approaches  have  been 
developed  for  this  problem,  including  black  box  and 
grey  box  techniques  (Bayard  and  Eslami  1984, 
Pajunen  1985,  Hunter  and  Korenberg  1986,  Korenberg 
and  Hunter  1986,  Hasiewicz  1987,  Chen  and  Fassois 
1992,  1997,  Greblicki  1992,  1994,  1997,  1998, 

Westwick  and  Kearney  1992,  Wigren  1994,  Westwick 
and  Verhaegen  1996,  Bai  1998,  Lovera  et  al.  2000, 
Van  Pelt  and  Bernstein  2000,  Lacy  and  Bernstein 
2001,  Lacy  et  cd.  2001,  Nelles  2001).  The  grey  box  case 
includes  the  identification  of  block-structured  models, 
such  as  the  Hammerstein  model  (linear  system  with 
input  non-linearity)  and  Weiner  model  (linear  system 
with  output  non-linearity)  (Haber  and  Keviczky 
1999  a,  b). 

This  paper  is  concerned  with  identifying  Wiener 
systems  under  more  general  assumptions  than  have 
been  previously  considered.  Many  methods  for  Wiener 
system  identification  require  the  non-linearity  to  be 
known,  invertible,  differentiable,  odd  or  require 
specially  designed  input  sequences.  In  particular,  the 
Wiener  identification  problem  has  been  considered  in 
Brillinger  (1970),  Pajunen  (1985),  Hasiewicz  (1987), 
Greblicki  (1992,  1994,  1997,  1998),  Westwick  and 
Kearney  (1992),  Wigren  (1994),  Westwick  and 
Verhaegen  (1996),  Bai  (1998),  and  Lovera  et  al.  (2000) 
under  the  assumption  that  the  non-linearity  is  unknown 
but  one-to-one.  This  assumption  simplifies  the  problem 
considerably  since  the  inverse  system  can  be  viewed  as  a 
Hammerstein  system  wherein  the  input  to  the  non- 
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linearity  is  measured.  If  the  non-linearity  is  known  but 
not  one-to-one,  then  identification  is  possible  by  first 
generating  a  candidate  set  of  signals  at  the  output  of 
the  linear  system  (Bayard  and  Eslami  1984,  Lacy  et  al. 
2001).  These  methods  are  applicable  even  if  the  output 
non-linearity  is  a  step  function  in  which  case  the  output 
assumes  at  most  two  distinct  values  (Lacy  et  al.  2001).  If 
the  input  sequence  can  be  chosen  freely,  the  frequency 
content  of  the  input  sequence  can  be  selected  such  that 
the  effect  of  the  non-linearity  can  be  derived  from  the 
frequency  content  of  the  output  (Pintelon  and 
Schoukens  2001). 

In  the  present  paper  we  consider  Wiener  system 
identification  in  which  the  output  non-linearity  is  both 
unknown  and  not  necessarily  one-to-one.  In  this  case 
the  goal  is  to  simultaneously  identify  both  the  linear 
system  dynamics  and  the  non-linearity  despite  the  non¬ 
in' vertibility  of  the  output  non-linearity.  To  do  this 
we  assume  that  the  non-linearity  can  be  represented  as 
a  finite  sum  of  polynomials.  We  use  multi-index 
notation  (Evans  1998,  Dunkl  and  Xu  2001)  to  expand 
this  polynomial  and  write  the  output  as  a  linear-in¬ 
parameters  sum  of  known  terms  with  unknown  coeffi¬ 
cients.  These  coefficients  consist  of  products  of  the 
system  parameters.  We  then  present  several  methods 
for  extracting  the  system  parameters.  This  approach  of 
expanding  the  polynomial  output  non-linearity  requires 
that  the  linear  system  output  depend  only  on  past 
inputs,  that  is,  the  linear  dynamics  are  assumed  to  be 
FIR. 

This  paper  is  organized  as  follows.  In  §2  we  define 
the  problem  and  list  the  assumptions.  In  §3  we  intro¬ 
duce  notation  used  throughout  the  paper.  In  §3.1  we 
present  an  algebraic  solution.  In  §3.2  we  present  a  sol¬ 
ution  based  on  a  singular  value  decomposition.  In  §3.3 
we  present  a  solution  based  on  a  multi-dimensional  sin¬ 
gular  value  decomposition  (Andersson  and  Bro  2000). 
In  §3.4  we  present  an  optimality  approach  based  on  a 
prediction-error  cost  function.  In  §4  we  apply  all  of 
these  methods  to  an  example  to  illustrate  their  imple¬ 
mentation  and  compare  their  effectiveness. 
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Figure  1.  Wiener  system. 


2.  Problem  description 

Here  we  study  the  identification  of  a  single-input, 
single-output,  linear  time-invariant  system  with  finite 
impulse  response  whose  output  is  measured  through  a 
static  non-linearity.  This  system,  which  is  represented  in 
figure  1,  is  modelled  by 


3.  Wiener  identification 

Using  (3),  we  rewrite  equation  (2)  as 

p  p  /  m 

z(k )  =  E  CiyW  =  E c'  E  hAk  -/> 

i= o  i=o  \y=o 


1=0 


M =i 


l«l=0 


=  4 i(k)Te 

where 


(4) 


y{k)  =  ( hQu(k )  +  h\u(k  ~  1)  +  •  •  •  +  hmu(k  —  m) 

m 

=  E  hiu(k  —  /)  (1) 

i= 0 

z(k)  =  Af{y(k))  (2) 

where  u  is  the  input  to  the  system,  y  is  the  unmeasured 
output  of  the  linear  system,  /?,  are  Markov  parameters, 
and  r  is  the  measured  output  of  the  non-linearity.  We 
assume  that  the  non-linear  function  A/" :  1F5  — 1F5  is  a 
polynomial  of  the  form 


Af{y{k))  =  Ec/jW'  (3) 

(=0 

If  A f  is  not  a  polynomial,  then  (3)  can  be  regarded  as  an 
approximation.  We  assume  that  the  order  m  of  the  FIR 
dynamics  and  the  degree  p  of  the  polynomial  A f  are 
known.  The  non-linearity  Af  is  otherwise  unknown 
and  not  necessarily  one-to-one.  In  the  practical  situation 
m  and  p  are  not  known.  Also,  it  may  be  difficult  to 
obtain  satisfactory  results  if  these  parameters  are  under¬ 
estimated.  However,  if  upper  bounds  on  these  par¬ 
ameters  are  known,  then  the  bounds  can  be  used  in  (1) 
and  (3)  at  the  expense  of  increasing  the  computational 
complexity.  Alternatively,  the  values  for  p  and  Af  can  be 
incremented  until  satisfactory  performance  is  achieved, 
at  the  expense  of  increasing  the  computational  load  at 
each  increment. 

The  identification  problem  is  to  estimate  the  coeffi¬ 
cients  hj  and  ct  using  l  measurements  of  u  and  z.  We 
adopt  a  two-stage  approach.  First,  we  solve  a  least 
squares  problem  to  obtain  an  estimate  0  of  a  vector  0 
whose  entries  are  the  unknown  parameters  and  products 
of  the  unknown  parameters.  Next,  we  present  several 
techniques  that  use  0  to  estimate  the  individual 
unknown  parameters.  In  addition,  we  minimize  a 
prediction  error  cost  function  to  further  refine  the 
parameter  estimates  and  compare  to  the  direct 
approaches. 


h  =  [ho  hi  ■■■  h,„}T  G  0r,+1  (5) 

a  =  [ax  a2  •••  am+1]T  G  Nq+1  (6) 

v(k)  =  [ u(k )  u(k  —  1)  •  •  •  u(k  —  m)\T  G  IK"'+1  (7) 

0=[C|a|A“]  |a|<„e^+1  (8) 


<t>{k)  = 


a! 


J  H</7 


(9) 


and  a  G  Nq,+1  is  a  multi-index  whose  order  is  m+  1, 
where  N0  is  the  set  of  positive  integers  and  zero. 
A  multi-index  is  a  vector  whose  components  are  non¬ 
negative  integers  (see  Evans  1998,  Dunkl  and  Xu  2001). 
We  define 


ra+1 

H  =  a\  +  a2  H - 1"  am+ 1  =  ^  ai  (10) 

i=l 


ra+1 

a\  =  (ai!)(o20  ' ' '  (am+iD  =  ]Ja!! 

1=1 


ra+1 

v(k)a  =  Vi  ( k)a'  v2(kr  ■  ■  ■  Vm+l  0 k)a =  n  Viiky  (12) 

1=1 


ra+1 


if  =  hfhf  ■  ■  ■  tiffl  =  n  K 


(13) 


i—  1 


The  number  of  multi-indices  of  order  m  +  1  of  fixed 
absolute  value  i  is  given  by 

(m  +  /)! 


^*ra+l  _ 


m  + 1 
i 


m  + 1 
m 


m\i\ 


(14) 


and  the  number  of  multi-indices  of  order  m  +  1  of  abso¬ 
lute  value  less  than  or  equal  to  p  is  given  by 


D 


fn+\  _  ^2  C'"+X  =  ( m  +  0 ■ 


7=0 


7=0 


m\i\ 


(15) 


We  need  to  define  an  order  relation  for  multi-indices  of 
the  same  order.  Let  a,  fj  G  N'0”+1  be  multi-indices.  If 
|a|  >  \/3\  then  a  >  (3.  If  a  and  (3  have  the  same  order, 
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|a|  =  \(j\,  then  we  choose  the  standard  dictionary  order¬ 
ing.  The  notation 

’  l/(a)]|a|=o" 

[/[(a)]|a|=l 

[f(a)]\a\<p=  .  (16) 

_  \f(a)\\a\=p  _ 

denotes  the  column  vector  whose  components  are / eval¬ 
uated  at  every  multi-index  a  such  that  |a|  <  p.  The  com¬ 
ponents  are  ordered  according  to  the  above  ordering 
scheme.  This  vector  has  Z>",+  1  components.  Thus  0  and 
(j>(k)  have  Z)”,+  l  components. 

To  estimate  0  we  rewrite  (4)  as 

x  =  <PTe  (17) 

where 

z  =  [z(m  +  1)  •••  :(<)]Tef"  (18) 

<Z>  =  +1)  •  •  •  </>(f)  G  uD"+'*t~m  (19) 

We  assume  <P<PT  is  non-singular,  which  is  a  persistency 
of  excitation  condition  that  requires  l  >  Z)”'+1  +  m. 
Then  we  calculate  the  least  squares  estimate  0  of  0 
given  by 

0=(<M>t)-‘$z  (20) 

Next  we  develop  several  methods  for  obtaining  esti¬ 
mates  c  and  h  based  on  0.  Note  that  an  arbitrary  scaling 
and  its  reciprocal  can  be  applied  to  the  linear  system  and 
the  output  non-linearity.  We  remove  this  ambiguity  by 
introducing  a  normalization  constraint,  thereby  select¬ 
ing  a  single  system  from  a  class  of  equivalent  systems. 
We  can  normalize  c  and  h  by  setting  c,-  =  a,  hj  =  a, 
||c||  =  a ,  ||/;||  =  a  or  various  other  constraints. 

3.1.  Direct  solve 

We  have  m  +  p  +  2  unknown  parameters  in  h  and  c, 
one  normalization  constraint  and  D”'+1  equations  in 
terms  of  0.  We  can  normalize  as  above,  then  choose 
m  +  p  +  1  independent  equations.  These  m  +  p  +  1 
equations  must  also  be  independent  of  the  constraint 
equation,  which  constitutes  equation  number 
m  +  p  +  2.  A  symbolic  manipulator  such  as 
Mathematica  can  be  used  to  invert  these  non-linear 
equations  and  obtain  estimates  h  and  c  of  h  and  c. 

3.2.  SVD 

To  begin,  c0  can  be  estimated  directly  by 

c0  =  0(1) 


To  estimate  the  remaining  components  of  c  and  h,  we 


arrange  the  components  of  0  into 
A(0)  G  Um+lxDp-' ,  where 

the  matrix 

A(0)  =  #T 

(22) 

and 

l/>=  [c\a\+lha]H<p 

(23) 

Then  we  calculate  the  singular  value  decomposition 

A(0)  =  USVT 

(24) 

to  obtain  the  estimates 

h  =  aS(  1,1)F(1,:) 

(25) 

a 

(26) 

where  the  scalar  a  ^  0  selects  the  normalization  con¬ 
straint.  Finally,  we  extract  the  non-linearity  coefficients 
dj  from  ip.  Specifically,  c t  is  given  directly  by  V’(l),  while 
the  remaining  coefficients  are  calculated  using  least 
squares  estimation  and  h. 

3.3.  Multi-dimensional  SVD 

First,  we  define  the  tensors  A0  G  IR,  Aj 

g  x-)=1or+1 

^o(0)  =  c0 

(27) 

Ai(0)  =  c\h  =  hAl 

(28) 

A2(0)  =  c2h  o  h  =  oj=1I, Ai 

(29) 

A}(0)  =  c2h  o  h  o  h  =  o-=i  hA 

(30) 

AP{0)  =  cp  of=1  h  =  o pi=lhAp 

(31) 

where 

hA„  =  c)!nh 

(32) 

We  use  a  multi-dimensional  singular  value  decomposi¬ 
tion  (Andersson  and  Bro  2000)  to  obtain  the  estimate 
hA.  of  hA..  To  do  this  we  note  that 

[hAl  >U2  ■  ■  ■  hAp }  =  hdT 

(33) 

where 

d=[c  i  cf  •••  c)!”]1 

(34) 

Flence  we  compute  the  singular  value  decomposition 

[hAi  ■  ■  ■  hAp\  =  USVT 

(35) 

(21) 

to  obtain  the  estimates 
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A  =  <rS(l,l)C7(:,l) 

(36) 

<9/pe 

d  =  -V(:,  1) 

(37) 

dcfcj 

a 

9/pe 

Cj  =  d\ 

(38) 

dcjdhj 

^eSe  S  e  Kwf 

\a\=i  |or|=/  ^  k=m+ 1 


~2i-  E  i  E  vW° 

|o|=i  fc=m+l 


(46) 


where  cr  f  0  selects  the  normalization  constraint. 

3.4.  Prediction  error  cost  function 

Consider  the  prediction  error  cost  function 

Jpfcjt)  =  ||z-z||  (39) 

where  z  is  the  response  of  the  estimated  system,  that  is 

P  /  m  \ 

z(k)  =  E E  hAk  ■■/)  (4°) 

!=0  \/=0  / 

Hence 

4(c,A)  =  (*(*)  -  E  )  (41) 

k— m+\  \  |a|— 0  J 


z(k)) 


m 

f3\ 


v(k)' 


(47) 


and 


c>4 

dhjdhj 


"2E^rcHa<  E  v(*0 

|o;|=0  *  /:=m+l 


(a,  -  6„)hr«-*(z(k)  -  z(/r))  -  /*“-«' 


^  |/?|! 


xE^w%,r 


<1=0 


(48) 


where 


The  derivatives  of  z(k)  are 


®  =  y  M  v(jfc)“A“  =  a  y  E  v(*r 

9c, ■  a!  '  rv! 


I«M 


l“M 


.  a! 


9z(A:)  4  lal'  ,-Q 

44  =  ,?^v(")  ci“'a,A 


|a|=0 


(42) 

(43) 


where  c,-  is  the  ith  column  of  Hence 


J  1  if  i  =  j 

9  |0  else 


(49) 


Using  the  above  expressions,  we  implement  a  gradi¬ 
ent-based  optimization  algorithm  to  minimize  /pe  and 
obtain  estimates  c  and  h.  Assuming  q  f  0,  we  normalize 
by  letting  q  =  1,  and  thus  remove  it  from  the  optimiza¬ 
tion  problem. 


°df=2±  (z(t)-W)4^(4 

UC‘  k—m+l  y  |a|=i  ) 

=  ~2il  EuE  *(*)“(*(*)  -  z(k))  (44) 

\a\=i  m+1 

and 

=  (^)-^))f-EPv(^)“cwa4-ei| 

k—m+ 1  y  |a|=0  a‘  / 

=  -2E^M^“~ei  E  <k)a{z{k)-z{k)) 

|q|=0  ‘  k—m+l 


4.  Example 

Let  m=  1,  p  =  3,  £  =  213  =  8192,  A  =  [2  1]T, 
c  =  [—20  1  10  1]T,  and  Jf{y)  =  -20  +y  +  lOy- 

y3,  which  is  shown  in  figure  2.  Thus 

y(£)  =  h0u(k)  4-  h\u{k  —  1)  (50) 

z(k)  =  c0  4-  C\h\ufk  —  1)  +  cftQu{k)  4-  c2h\u{k  —  l)2 

T  ’2.c2hohiu(k')u(k  —  1)4-  c2 u{!f 

4“  c2h\u(k  —  1)  4-  ’ic2hftQu{k)u{k  —  1) 

4“  3cih0hiu(k)  u(k  —  1)4-  c^h^u^kf  4-  vv {k) 


(45)  =  f{k)T0  +  w(k) 


(51) 


The  second  derivatives  are  given  by 


where 
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JVC!/) 


y(*: 

Figure  2.  A f(y). 


0  =  \c\a\ha\\a\<p 

—  [co  f  i  //  c,/i0  C2/1?  Cihfji  Cihl  c}h]  c3h0h\  C^h^h  c^lil] 

(52) 

o  =  fM 

L  «!  J  |a|<J7 

=  [1  w(fc  —  1)  u(k)  u(k  —  l)2  2u(k)u(k  —  1)  m(A)2  u(k  —  l)3 

3u(k)u(k  —  l)2  3u(k)2u(k  —  1)  /*(£)3]T  (53) 


and  u’(/r)  is  a  realization  of  a  zero-mean  Gaussian  pro¬ 
cess  such  that  the  signal  to  noise  ratio 


s/w= =  10 
HI 

(54) 

Then  (17)  is 

replaced  by 

z  =  (pT0  +  W 

(55) 

where 

z  =  [z(  2) 

•••  mT 

(56) 

0  =  [cp(  2) 

...  m 

(57) 

w  =  [w(m 

+  1)  •••  w(t)\T  =  [w( 2)  •• 

OO 

H 

T 

We  assume  is  non-singular,  and  estimate  the  par¬ 
ameter  vector  using  (20). 


4.1.  Direct  solve 

We  have  m  +  p  +  2  =  6  unknown  parameters,  one 
normalization  constraint,  and  Z>”,+1  =  10  equations. 
We  choose  to  normalize  by  setting  c\  =  1,  and  solve 
m+p  +  1  =  5  equations.  Using  0(1),  0(2),  0(3),  0(4), 
0(7)  and  ignoring  the  rest  of  0  we  obtain 


C0  =  0(1),  Cl  =  1,  £2  =  ^2  (59) 

0(2) 

G=S,  *0  =  0(2),  *i=0(3)  (60) 

0(2) 

Figure  3(a)  shows  h  for  100  simulations,  each  having  a 
different  realization  of  the  input  and  noise  sequences. 
Figure  3(b)  shows  the  identified  non-linearity. 

4.2.  SVD 

Flere  we  arrange  the  components  of  0  into  a  matrix 
that  is  also  an  outer  product  of  h  and  1)2.  Then  we  com¬ 
pute  the  singular  value  decomposition  of  this  matrix  to 
find  h  and  ip.  Finally,  we  extract  c  from  1 p.  First,  c0  can 
be  estimated  directly  as  in  (21).  To  find  the  remaining 
parameters,  we  arrange  the  components  of  0  as 

m  =  #T  =  hai+i*“]H</T 


Cl 

-  0(3) 

0(2)' 

c2h\ 

0(5) 

0(4) 

c2h0 

[ho  hi]  = 

0(6) 

0(5) 

c3*i 

0(8) 

0(7) 

C3*0*l 

0(9) 

0(8) 

1 

(NO 

-s: 

- 1 

.0(10) 

0(9) . 

We  calculate  the  singular  value  decomposition  of  A  (0) 
as  in  (24).  Then  we  obtain  h  and  ip  from  (25)  and  (26). 
Next,  we  extract  the  non-linearity  coefficients  c,-  from  ip. 
Ci  is  given  directly  by  ip(l),  and  we  calculate  the  remain¬ 
ing  Cj  using  least  squares  estimation.  In  this  case  the 
normalization  constraint  is  selected  by  the  choice  of 
the  scalar  a.  We  choose  to  normalize  by  setting 
<r=  \J(  1, 1)  such  that  c\  =  ip(  1)  =  1.  In  figure  3(c)  we 
plot  h  for  100  simulations.  In  figure  3(d)  we  plot  the 
identified  non-linearity. 


4.3.  Multi-dimensional  SVD 

We  arrange  the  elements  of  0  into  several  matrices  as 
in  (27)  (3 1 ),  corresponding  to  c,-  and  the  z'th  power  of  h 


Aq  —  eg  —  0(1) 


A-i  — 

J^2  = 


C\h  =  C[ 


'V 

-0(3)- 

.hi. 

.0(2). 

'hQ' 

c2 

[ho  h\ 

hi. 

=  hA 


0(6)  0(5) 
0(5)  0(4) 


=  hA  o  hA 


(62) 

(63) 


(64) 
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(a)  Direct  Solve  Solutions:  h\  vs  ho  (b)  Direct  Solve  Solutions:  If 


(c)  SVD  Solutions:  h\  vs  ho  (d)  SVD  Solutions:  Jf 


(e)  Multi-Dimensional  SVD  Solutions:  (f)  Multi-Dimsional  SVD  Solutions:  Jf 
hi  vs  ho 


(g)  Optimization  Solutions:  hi  vs  ho  (h)  Optimization  Solutions:  Jf 
Figure  3.  Simulation  results. 
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A  3  =  c3/t  0/10/1  =  hA}  o  hA}  o  hA} 


’  h30 

hoh\ 

>10) 

0(9)' 

h\h\ 

h0hi 

= 

0(9) 

0(8) 

Iiq/ii 

h0hi 

0(9) 

0(8) 

.  hah\ 

h\  . 

.  m 

0(  7). 

We  use  a  multi-dimensional  singular  value  decomposi¬ 
tion  (Andersson  and  Bro  2000)  to  estimate  the  scaled 
Markov  vector  hAj.  Next,  we  arrange  these  results  in  a 
matrix  and  compute  the  singular  value  decomposition  as 
in  (35)  to  estimate  h ,  d ,  and  c,  see  (36)  (38).  In  this  case 
the  normalization  constraint  is  enforced  by  how  we 
choose  a.  We  choose  to  normalize  by  setting 
a=  VA(l,  1)  resulting  in  c {  =  1.  In  figure  3(e)  we  plot 
h  for  100  simulations.  In  figure  3(f)  we  plot  the  identified 
non-linearity. 

4.4.  Prediction  error  cost  function 

We  minimize  the  standard  prediction  error  cost  func¬ 
tion,  /pe(c,  h)  =  ||z  —  z||,  initialized  with  one  of  the  three 
algorithms  discussed  earlier.  Using  the  derivatives  given 
previously  we  obtain  c*  and  h* .  Although  the  results  of 
the  optimization  are  generally  independent  of  the 
method  used  to  obtain  the  initial  estimate,  the  number 
of  iterations  needed  by  the  optimization  routine  to 
converge  usually  depends  on  the  accuracy  of  the  initial 
estimate.  We  used  lsqnonlin  in  the  Matlab  optimi¬ 
zation  toolbox  to  minimize  the  function.  In  figure  3(g) 
we  plot  h  for  100  simulations.  In  figure  3(h)  we  plot  the 
identified  non-linearity.  Due  to  noise  effects,  the  value  of 
the  cost  function  evaluated  at  the  optimal  parameters  is 
generally  less  than  the  cost  function  evaluated  at  the  true 
parameters  Jpe(c*,h*)  <  Jpe(c, h). 

4.5.  Discussion 

The  four  methods  presented,  namely  direct  algebraic 
solution,  singular  value  decomposition,  multi-dimen¬ 
sional  singular  value  decomposition  and  prediction- 
error  minimization,  produced  solutions  of  varying  accu¬ 
racy.  Solving  some  of  the  equations  for  the  unknown 
parameters  generally  produced  poor  quality  estimates, 
as  compared  to  the  other  three.  Since  the  direct  algebraic 
solution  method  uses  only  a  fraction  of  the  entries  of  0 
to  compute  the  estimates,  it  is  less  robust  than  the  other 
three  methods.  In  addition,  the  direct  algebraic  solution 
method  requires  user  interaction  to  select  which  equa¬ 
tions  to  solve. 

The  singular  value  decomposition  approch  produced 
estimates  that  were  close  to  the  ones  obtained  by  the 
prediction  error  minimization,  but  relied  only  on  a  sin¬ 


gular  value  decomposition  and  some  subsequent  least 
squares  steps.  This  method  is  the  simplest  to  implement 
both  for  the  user  and  numerically. 

The  multi-dimensional  singular  value  decomposition 
approach  produced  estimates  comparable  to  the  first 
singular  value  decomposition  approach,  but  the  multi¬ 
dimensional  singular  value  decomposition  is  more  com¬ 
plex  to  implement  than  the  two-dimensional  singular 
value  decomposition  approach. 

The  prediction-error  minimization  approach  is  the 
most  complex  to  implement.  For  best  results,  it  should 
be  initialized  using  one  of  the  previous  methods. 
Flowever,  it  produced  the  best  estimates  of  both  the 
linear  dynamics  and  the  output  non-linearity. 


5.  Conclusion 

This  paper  presented  four  methods  for  identifying 
FIR  Wiener  systems  with  polynomial  non-linearities. 
We  presented  three  methods  for  simultaneous  direct 
estimation  of  the  non-linearity  and  linear  dynamics, 
and  a  prediction  error  optimization  method. 

Many  authors  have  studied  Wiener  system  identifi¬ 
cation  under  the  assumption  that  the  non-linearity  is 
unknown  but  one-to-one  (Brillinger  1970,  Pajunen 
1985,  Hasiewicz  1987,  Greblicki  1992,  1994,  1997, 
Westwick  and  Kearney  1992,  Wigren  1994,  Westwick 
and  Verhaegen  1996,  Bai  1998,  Lovera  et  al.  2000). 
Other  methods  for  Wiener  system  identification  require 
the  non-linearity  to  be  known,  invertible,  monotonic, 
odd,  even,  or  require  the  use  of  specially  designed 
input  sequences.  In  this  paper  we  require  the  non-lin¬ 
earity  to  be  polynomial  and  the  linear  dynamics  to  have 
finite  impulse  response.  These  two  assumptions  are  prac¬ 
tical  in  that  many  non-linearities  can  be  approximated 
with  polynomials,  and  that  many  systems  with  infinite 
impulse  response  can  be  approximated  with  finite 
impulse  response  dynamics. 

Future  work  will  focus  on  extending  the  method  in 
three  directions:  first,  the  identification  of  sandwich  non¬ 
linear  systems,  i.e.  systems  with  both  input  and  output 
non-linearities;  second,  the  identification  of  HR  Wiener 
systems;  and  finally,  identification  of  Wiener  systems 
with  non-polynomial  output  non-linearities.  While  the 
identification  of  sandwich  non-linear  systems  using  this 
approach  seems  tractable,  overcoming  the  FIR  and 
polynomial  assumptions  on  the  linear  dynamics  and 
output  non-linearity  appears  to  be  more  challenging. 
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